1 Introduction

2 Preparation

source("Port_opt_semi-G-structure-functions.R")
## 
## Attaching package: 'CVXR'
## The following object is masked from 'package:matrixcalc':
## 
##     vec
## The following object is masked from 'package:stats':
## 
##     power
## Found more than one class "atomicVector" in cache; using the first, from namespace 'Matrix'
## Also defined by 'Rmpfr'
## Found more than one class "atomicVector" in cache; using the first, from namespace 'Matrix'
## Also defined by 'Rmpfr'
## Found more than one class "atomicVector" in cache; using the first, from namespace 'Matrix'
## Also defined by 'Rmpfr'
## Found more than one class "atomicVector" in cache; using the first, from namespace 'Matrix'
## Also defined by 'Rmpfr'
## Found more than one class "atomicVector" in cache; using the first, from namespace 'Matrix'
## Also defined by 'Rmpfr'
## Found more than one class "atomicVector" in cache; using the first, from namespace 'Matrix'
## Also defined by 'Rmpfr'
## Found more than one class "atomicVector" in cache; using the first, from namespace 'Matrix'
## Also defined by 'Rmpfr'
## Found more than one class "atomicVector" in cache; using the first, from namespace 'Matrix'
## Also defined by 'Rmpfr'
## Found more than one class "atomicVector" in cache; using the first, from namespace 'Matrix'
## Also defined by 'Rmpfr'
## Found more than one class "atomicVector" in cache; using the first, from namespace 'Matrix'
## Also defined by 'Rmpfr'
## Found more than one class "atomicVector" in cache; using the first, from namespace 'Matrix'
## Also defined by 'Rmpfr'
## Found more than one class "atomicVector" in cache; using the first, from namespace 'Matrix'
## Also defined by 'Rmpfr'
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.5     ✓ purrr   0.3.4
## ✓ tibble  3.1.6     ✓ dplyr   1.0.7
## ✓ tidyr   1.1.4     ✓ stringr 1.4.0
## ✓ readr   2.1.1     ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter()    masks stats::filter()
## x dplyr::id()        masks CVXR::id()
## x purrr::is_vector() masks CVXR::is_vector()
## x dplyr::lag()       masks stats::lag()
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
sd.mat0 <- matrix(c(1,3, 
                    2,4), 2, 2, byrow = TRUE)
mu.vec <- matrix(c(1, 2), 2, 1)
rho.vec <- c(0.2,0.4)

3 Implementation

3.1 P: min variance

(A1 <- A.f(sd.mat0, rho.vec, par.weight = 1/2, i=1))
##      [,1] [,2]
## [1,]  5.0  2.6
## [2,]  2.6 10.0
(A2 <- A.f(sd.mat0, rho.vec, par.weight = 3/4, i=2))
##      [,1] [,2]
## [1,] 1.75 0.85
## [2,] 0.85 4.50
is.positive.definite(A1)
## [1] TRUE
is.positive.definite(A2)
## [1] TRUE
# Variables minimized over
x <- Variable(2)

# Problem definition
objective <- Minimize(quad_form(x, A2/2))
constraints <- list(x >= 0, sum(x) == 1)
#constraints <- list(sum(x) == 1)
prob <- Problem(objective, constraints)

# Problem solution
solution <- solve(prob)
solution$status
## [1] "optimal"
solution$value
## [1] 0.785989
solution$getValue(x)
##           [,1]
## [1,] 0.8021978
## [2,] 0.1978022
#average allocation
#achieve smallest variance uncertainty
#consider another problem
# Variables minimized over
x <- Variable(2)
t <- 1
# Problem definition
objective <- Minimize(quad_form(x, A1/2) - t*t(mu.vec)%*%x)
constraints <- list(x >= 0, sum(x) == 1)
#constraints <- list(sum(x) == 1)
prob <- Problem(objective, constraints)

# Problem solution
solution <- solve(prob)
solution$status
## [1] "optimal"
solution$value
## [1] 0.9102041
solution$getValue(x)
##           [,1]
## [1,] 0.6530612
## [2,] 0.3469388
(H.opt <- solution$getValue(quad_form(x, A1)))
## [1] 4.514286
(mu.opt <- solution$getValue(t(mu.vec)%*%x))
## [1] 1.346939

3.2 Q: efficient frontier

3.2.1 A1

#draw the efficient frontier 
t.seq <- seq(-5,5,.1)
H.opt.seq <- mu.opt.seq <- numeric(length(t.seq))
# x.opt.mat <- matrix(NA, 
#                     nrow = length(t.seq), ncol=2)
for (i in seq_along(t.seq)){
  t <- t.seq[i]
  x <- Variable(2)
# Problem definition
objective <- Minimize(quad_form(x, A1/2) - t*t(mu.vec)%*%x)
constraints <- list(x >= 0, sum(x) == 1)
#constraints <- list(sum(x) == 1)
prob <- Problem(objective, constraints)

# Problem solution
solution <- solve(prob)
#x.opt.mat[i,] <- as.numeric(solution$getValue(x))
H.opt.seq[i] <- solution$getValue(quad_form(x, A1))
mu.opt.seq[i] <- solution$getValue(t(mu.vec)%*%x)
}
#efficient frontier
plot(H.opt.seq, mu.opt.seq, type = "l", 
     main = "A1 opt problem")

3.2.2 A2

#draw the efficient frontier 
t.seq <- seq(-5,5,.1)
H.opt.seq <- mu.opt.seq <- numeric(length(t.seq))
for (i in seq_along(t.seq)){
  t <- t.seq[i]
  x <- Variable(2)
# Problem definition
objective <- Minimize(quad_form(x, A2/2) - t*t(mu.vec)%*%x)
constraints <- list(x >= 0, sum(x) == 1)
#constraints <- list(sum(x) == 1)
prob <- Problem(objective, constraints)

# Problem solution
solution <- solve(prob)
H.opt.seq[i] <- solution$getValue(quad_form(x, A1))
mu.opt.seq[i] <- solution$getValue(t(mu.vec)%*%x)
}
#efficient frontier
plot(H.opt.seq, mu.opt.seq, type = "l", 
     main = "A2 opt problem")

EF does depend on the measure we choose and also the tuning parameter (which depends on the preference of the user on the variance uncertainty).

3.3 EF: change par

3.3.1 A1

t.seq <- seq(-10, 10, .2)
a.seq <- seq(0.1, 0.9,.1)
H.opt.mat <- mu.opt.mat <- matrix(NA, 
                                  ncol = length(a.seq),
                                  nrow = length(t.seq))


#we can run parallel computing for this part
for (j in seq_along(a.seq)){
  a <- a.seq[j]
  A1 <- A.f(sd.mat0, rho.vec, par.weight = a, i=1)
  cat("j = ",j, "\n")
for (i in seq_along(t.seq)){
  t <- t.seq[i]
  x <- Variable(2)
  # Problem definition
  objective <- Minimize(quad_form(x, A1/2) - t*t(mu.vec)%*%x)
  constraints <- list(x >= 0, sum(x) == 1)
  #constraints <- list(sum(x) == 1)
  prob <- Problem(objective, constraints)

  # Problem solution
  solution <- solve(prob)
  H.opt.mat[i, j] <- solution$getValue(quad_form(x, A1))
  mu.opt.mat[i, j] <- solution$getValue(t(mu.vec)%*%x)
}

}

#plot the matrix
#get the largest range
H.lim <- c(min(H.opt.mat), max(H.opt.mat))
mu.lim <- c(min(mu.opt.mat), max(mu.opt.mat))

for (j in seq_along(a.seq)){
  if (j==1){
    plot(H.opt.mat[,j],mu.opt.mat[,j], col=j,lty=j,
         type = "l", xlab = "H.opt", ylab = "mu.opt", 
         xlim = H.lim, ylim = mu.lim, 
         main = "Changing tuning par for P2.1")
  } else {
    lines(H.opt.mat[,j],mu.opt.mat[,j], col=j,lty=j,
          type = "l")
  }
}

3.3.2 A2

t.seq <- seq(-10, 10, .2)
a.seq <- seq(2/3, 0.99, length.out = 9)
H.opt.mat2 <- mu.opt.mat2 <- matrix(NA, 
                                  ncol = length(a.seq),
                                  nrow = length(t.seq))
#we can run parallel computing for this part
#seq_along(a.seq)
for (j in seq_along(a.seq)){
  a <- a.seq[j]
  A1 <- A.f(sd.mat0, rho.vec, par.weight = a, i=1)
  cat("j = ",j, "\n")
for (i in seq_along(t.seq)){
  t <- t.seq[i]
  x <- Variable(2)
  # Problem definition
  objective <- Minimize(quad_form(x, A1/2) - t*t(mu.vec)%*%x)
  constraints <- list(x >= 0, sum(x) == 1)
  #constraints <- list(sum(x) == 1)
  prob <- Problem(objective, constraints)

  # Problem solution
  solution <- solve(prob)
  H.opt.mat2[i, j] <- solution$getValue(quad_form(x, A1))
  mu.opt.mat2[i, j] <- solution$getValue(t(mu.vec)%*%x)
}

}

#plot the matrix
#get the largest range
H.lim <- c(min(H.opt.mat2), max(H.opt.mat2))
mu.lim <- c(min(mu.opt.mat2), max(mu.opt.mat2))

for (j in seq_along(a.seq)){
  if (j==1){
    plot(H.opt.mat2[,j],mu.opt.mat2[,j], col=j,lty=j,
         type = "l", xlab = "H.opt", ylab = "mu.opt", 
         xlim = H.lim, ylim = mu.lim, 
         main = "Changing tuning par for P2.2")
  } else {
    lines(H.opt.mat2[,j],mu.opt.mat2[,j], col=j,lty=j,
          type = "l")
  }
}

4 Extended cases

4.1 Objective function

#consider the rho as negative
#consider this problem from theoretical level
#check whether the optimal weight depend on the tuning parameter
#plot the range of the functions 
#H is the objective function
#draw a 3D surface plot
sig1 <- seq(0.1, 4, .01)
sig2 <- seq(2, 5, .01)

z <- outer(sig1, sig2, H)

contour(sig1, sig2, z, xlab="sig1", ylab="sig2", xlim = range(sig1), ylim = range(sig2))

#abline(a=0, b=0.5, col=2, xlim = range(sig1))
#abline(a=0, b=1/0.5, col=2, ylim = range(sig2))
which.min(z)
## [1] 16
which.max(z)
## [1] 391
plot(function(x) H(x, sig2=2), from = min(sig1), to = max(sig1), 
     xlab = "sig1")

plot(function(x) H(sig1=2, sig2=x), from = min(sig2), to = max(sig2), 
     xlab = "sig2")

x <- sig1
y <- sig2
fig <- plot_ly(
  type = 'surface',
  contours = list(
    z = list(show = TRUE, start = 0, end = max(z), size = 0.04)),
  x = ~x,
  y = ~y,
  z = ~z)
fig <- fig %>% layout(
    scene = list(
      xaxis = list(nticks = 20),
      zaxis = list(nticks = 20),
      camera = list(eye = list(x = 0, y = -1, z = 0.5)),
      aspectratio = list(x = .9, y = .8, z = .6)))
fig
#our current guess is, for 
#check different cases
#the objective function may no longer be a convex one
#use range.set
range.set(x1=0.2, print.ind = TRUE, rho.intvl = c(-0.8, -0.4))
## i =  1 , sig1 =  0.5 , sig2 =  1 , rho =  -0.8 , H =  0.522 
## i =  2 , sig1 =  0.5 , sig2 =  2 , rho =  -0.8 , H =  2.314 
## i =  3 , sig1 =  1 , sig2 =  1 , rho =  -0.8 , H =  0.424 
## i =  4 , sig1 =  1 , sig2 =  2 , rho =  -0.8 , H =  2.088 
## i =  5 , sig1 =  0.5 , sig2 =  1 , rho =  -0.4 , H =  0.586 
## i =  6 , sig1 =  0.5 , sig2 =  2 , rho =  -0.4 , H =  2.442 
## i =  7 , sig1 =  1 , sig2 =  1 , rho =  -0.4 , H =  0.552 
## i =  8 , sig1 =  1 , sig2 =  2 , rho =  -0.4 , H =  2.344
## [1] 0.424 2.442
range.set(x1=0.5, print.ind = TRUE, rho.intvl = c(-0.8, -0.4))
## i =  1 , sig1 =  0.5 , sig2 =  1 , rho =  -0.8 , H =  0.1125 
## i =  2 , sig1 =  0.5 , sig2 =  2 , rho =  -0.8 , H =  0.6625 
## i =  3 , sig1 =  1 , sig2 =  1 , rho =  -0.8 , H =  0.1 
## i =  4 , sig1 =  1 , sig2 =  2 , rho =  -0.8 , H =  0.45 
## i =  5 , sig1 =  0.8 , sig2 =  1 , rho =  -0.8 , H =  0.09 
## i =  6 , sig1 =  0.8 , sig2 =  2 , rho =  -0.8 , H =  0.52 
## i =  7 , sig1 =  0.5 , sig2 =  1 , rho =  -0.4 , H =  0.2125 
## i =  8 , sig1 =  0.5 , sig2 =  2 , rho =  -0.4 , H =  0.8625 
## i =  9 , sig1 =  1 , sig2 =  1 , rho =  -0.4 , H =  0.3 
## i =  10 , sig1 =  1 , sig2 =  2 , rho =  -0.4 , H =  0.85 
## i =  11 , sig1 =  0.8 , sig2 =  1 , rho =  -0.4 , H =  0.25 
## i =  12 , sig1 =  0.8 , sig2 =  2 , rho =  -0.4 , H =  0.84 
## i =  13 , sig1 =  0.8 , sig2 =  1 , rho =  -0.4 , H =  0.25 
## i =  14 , sig1 =  0.8 , sig2 =  2 , rho =  -0.4 , H =  0.84
## [1] 0.0900 0.8625
range.set(x1=0.8, print.ind = TRUE, rho.intvl = c(-0.8, -0.4))
## i =  1 , sig1 =  0.5 , sig2 =  1 , rho =  -0.8 , H =  0.072 
## i =  2 , sig1 =  0.5 , sig2 =  2 , rho =  -0.8 , H =  0.064 
## i =  3 , sig1 =  0.5 , sig2 =  1.6 , rho =  -0.8 , H =  0.0576 
## i =  4 , sig1 =  1 , sig2 =  1 , rho =  -0.8 , H =  0.424 
## i =  5 , sig1 =  1 , sig2 =  2 , rho =  -0.8 , H =  0.288 
## i =  6 , sig1 =  1 , sig2 =  1.6 , rho =  -0.8 , H =  0.3328 
## i =  7 , sig1 =  0.5 , sig2 =  1 , rho =  -0.4 , H =  0.136 
## i =  8 , sig1 =  0.5 , sig2 =  2 , rho =  -0.4 , H =  0.192 
## i =  9 , sig1 =  0.5 , sig2 =  1.6 , rho =  -0.4 , H =  0.16 
## i =  10 , sig1 =  0.5 , sig2 =  1.6 , rho =  -0.4 , H =  0.16 
## i =  11 , sig1 =  1 , sig2 =  1 , rho =  -0.4 , H =  0.552 
## i =  12 , sig1 =  1 , sig2 =  2 , rho =  -0.4 , H =  0.544 
## i =  13 , sig1 =  1 , sig2 =  1.6 , rho =  -0.4 , H =  0.5376 
## i =  14 , sig1 =  1 , sig2 =  1.6 , rho =  -0.4 , H =  0.5376
## [1] 0.0576 0.5520
#plot the maximum variance
x1.seq <- seq(-1, 2, .01)
#x1.seq <- seq(0, 1, .01)
rho.intvl1 <- c(-0.8,-0.4)
h.seq <- sapply(x1.seq, function(x){
  re <- range.set(x1=x, rho.intvl = rho.intvl1)
  re[2]
})
plot(x1.seq, h.seq, type = "l", ylab = "max var", 
     main = paste0("rho.intvl = [", rho.intvl1[1], 
                   ",", rho.intvl1[2], "]"))

#the max variance
#it is still a convex function
#not only for positive x1
rho.intvl1 <- c(-0.8,0.4)
h.seq <- sapply(x1.seq, function(x){
  re <- range.set(x1=x, rho.intvl = rho.intvl1)
  re[2]
})
plot(x1.seq, h.seq, type = "l", ylab = "max var", 
     main = paste0("rho.intvl = [", rho.intvl1[1], 
                   ",", rho.intvl1[2], "]"))

rho.intvl1 <- c(-0.8,0.4)
h.seq <- sapply(x1.seq, function(x){
  re <- range.set(x1=x, rho.intvl = rho.intvl1)
  re[1]
})
plot(x1.seq, h.seq, type = "l", ylab = "min var", 
     main = paste0("rho.intvl = [", rho.intvl1[1], 
                   ",", rho.intvl1[2], "]"))

rho.intvl1 <- c(-0.2,0.2)
h.seq <- sapply(x1.seq, function(x){
  re <- range.set(x1=x, rho.intvl = rho.intvl1)
  re[2]
})
plot(x1.seq, h.seq, type = "l", ylab = "max var", 
     main = paste0("rho.intvl = [", rho.intvl1[1], 
                   ",", rho.intvl1[2], "]"))

#plot the min variance
x1.seq <- seq(-1, 2, .01)
#x1.seq <- seq(0, 1, .01)
rho.intvl1 <- c(-0.8,-0.4)
h.seq <- sapply(x1.seq, function(x){
  re <- range.set(x1=x, rho.intvl = rho.intvl1)
  re[1]
})
plot(x1.seq, h.seq, type = "l", ylab = "min var", 
     main = paste0("rho.intvl = [", rho.intvl1[1], 
                   ",", rho.intvl1[2], "]"))

#the min variance
#it is not a convex function
#not only for positive x1
x1.seq <- seq(0.9, 1.1, .01)
rho.intvl1 <- c(-0.8,-0.4)
h.seq <- sapply(x1.seq, function(x){
  re <- range.set(x1=x, rho.intvl = rho.intvl1)
  re[1]
})
plot(x1.seq, h.seq, type = "l", ylab = "min var", 
     main = paste0("rho.intvl = [", rho.intvl1[1], 
                   ",", rho.intvl1[2], "]"))

x1.seq <- seq(-1, 2, .01)
rho.intvl1 <- c(-0.8, -0.5)
h.seq <- sapply(x1.seq, function(x){
  re <- range.set(x1=x, rho.intvl = rho.intvl1)
  re[2]-re[1]
})
plot(x1.seq, h.seq, type = "l", ylab = "max var", 
     main = paste0("rho.intvl = [", rho.intvl1[1], 
                   ",", rho.intvl1[2], "]"))

#It may not be convex
#but it should still have a global minimum

4.2 P: min variance

4.2.1 A1

#run the numerical optimization for this function
#so far we do not have a closed form solution

#direct minimization of H w.r.t x1

#both report the variance and current mean
E1 <- function(x, lam = 0.5, rho.intvl = c(-0.8, -0.5)){
  re <- range.set(x1=x, rho.intvl = rho.intvl)
  lam*re[2]+(1-lam)*re[1]
}

(opt.re <- optimize(E1, interval = c(-50,50)))
## $minimum
## [1] 0.7071823
## 
## $objective
## [1] 0.2370166
x1.opt <- opt.re$minimum

(x.opt <- c(x1.opt, 1-x1.opt))
## [1] 0.7071823 0.2928177
(obj.opt <- opt.re$objective)
## [1] 0.2370166

4.2.2 A2

E2 <- function(x, kap = 0.5, rho.intvl = c(-0.8, -0.5)){
  re <- range.set(x1=x, rho.intvl = rho.intvl)
  cen <- (re[2]+re[1])/2
  amb <- re[2]-re[1]
  kap*cen+(1-kap)*amb
}

(opt.re <- optimize(E2, c(-50, 50)))
## $minimum
## [1] 0.7173543
## 
## $objective
## [1] 0.3098996
x1.opt <- opt.re$minimum

(x.opt <- c(x1.opt, 1-x1.opt))
## [1] 0.7173543 0.2826457
(obj.opt <- opt.re$objective)
## [1] 0.3098996

4.3 Q: efficient frontier

  • change the objective function
  • subtract by \(tx\mu\)
  • then also draw the EF

4.3.1 A1

#draw the efficient frontier 
t.seq <- seq(-5,5,.1)
H.opt.seq <- mu.opt.seq <- numeric(length(t.seq))
mu.vec <- as.numeric(mu.vec)
# x.opt.mat <- matrix(NA, 
#                     nrow = length(t.seq), ncol=2)

rho.intvl1 <- c(-0.8, -0.5)
#objective function
E1 <- function(x, lam = 0.5, rho.intvl = rho.intvl1){
  re <- range.set(x1=x, rho.intvl = rho.intvl)
  lam*re[2]+(1-lam)*re[1]
}

for (i in seq_along(t.seq)){
  t <- t.seq[i]
  
  obj <- function(x) E1(x) - t*sum(mu.vec*x)
  opt.re <- optimize(obj, interval = c(-50,50))
  
  x1.opt <- opt.re$minimum
  x.opt <- c(x1.opt, 1-x1.opt)
  H.opt.seq[i] <- E1(x1.opt)
  mu.opt.seq[i] <- sum(mu.vec*x.opt)
}


#efficient frontier
plot(H.opt.seq, mu.opt.seq, type = "l", 
     main = paste0("Q1 opt problem with ", "rho.intvl = [", rho.intvl1[1],
                   ",", rho.intvl1[2], "]"))

4.3.2 A2

#draw the efficient frontier 
t.seq <- seq(-5,5,.1)
H.opt.seq <- mu.opt.seq <- numeric(length(t.seq))
mu.vec <- as.numeric(mu.vec)
# x.opt.mat <- matrix(NA, 
#                     nrow = length(t.seq), ncol=2)

#objective function
rho.intvl1 <- c(-0.8, -0.5)

E2 <- function(x, kap = 0.5, rho.intvl = rho.intvl1){
  re <- range.set(x1=x, rho.intvl = rho.intvl)
  cen <- (re[2]+re[1])/2
  amb <- re[2]-re[1]
  kap*cen+(1-kap)*amb
}

for (i in seq_along(t.seq)){
  t <- t.seq[i]
  
  obj <- function(x) E2(x) - t*sum(mu.vec*x)
  opt.re <- optimize(obj, interval = c(-50,50))
  
  x1.opt <- opt.re$minimum
  x.opt <- c(x1.opt, 1-x1.opt)
  H.opt.seq[i] <- E2(x1.opt) #to change 
  mu.opt.seq[i] <- sum(mu.vec*x.opt)
}


#efficient frontier
plot(H.opt.seq, mu.opt.seq, type = "l", 
     main = paste0("Q1 opt problem with ", "rho.intvl = [", rho.intvl1[1],
                   ",", rho.intvl1[2], "]"))

4.4 EF: change par

t.seq <- seq(-5, 5, .1)
a.seq <- seq(0.1, 0.9,.1)

4.4.1 A1

H.opt.mat <- mu.opt.mat <- matrix(NA, 
                                  ncol = length(a.seq),
                                  nrow = length(t.seq))

#we can run parallel computing for this part
for (j in seq_along(a.seq)){
  a <- a.seq[j]
  E1 <- function(x, lam = a, rho.intvl = c(-0.8, 0.5)){
  re <- range.set(x1=x, rho.intvl = rho.intvl)
  lam*re[2]+(1-lam)*re[1]
}
  #A1 <- A.f(sd.mat0, rho.vec, par.weight = a, i=1)
  cat("j = ",j, "\n")
for (i in seq_along(t.seq)){
  t <- t.seq[i]
  obj <- function(x) E1(x) - t*sum(mu.vec*x)
  opt.re <- optimize(obj, interval = c(-50,50))
  x1.opt <- opt.re$minimum
  x.opt <- c(x1.opt, 1-x1.opt)
  H.opt.mat[i, j] <- E1(x1.opt)
  mu.opt.mat[i, j] <- sum(mu.vec*x.opt)
}
}
## j =  1 
## j =  2 
## j =  3 
## j =  4 
## j =  5 
## j =  6 
## j =  7 
## j =  8 
## j =  9
#plot the matrix
#get the largest range
H.lim <- c(min(H.opt.mat), max(H.opt.mat))
mu.lim <- c(min(mu.opt.mat), max(mu.opt.mat))

for (j in seq_along(a.seq)){
  if (j==1){
    plot(H.opt.mat[,j],mu.opt.mat[,j], col=j,lty=j,
         type = "l", xlab = "H.opt", ylab = "mu.opt", 
         xlim = H.lim, ylim = mu.lim, 
         main = "Changing tuning par for Q1")
  } else {
    lines(H.opt.mat[,j],mu.opt.mat[,j], col=j,lty=j,
          type = "l")
  }
}

4.4.2 A2

H.opt.mat <- mu.opt.mat <- matrix(NA, 
                                  ncol = length(a.seq),
                                  nrow = length(t.seq))

#we can run parallel computing for this part
for (j in seq_along(a.seq)){
  a <- a.seq[j]
  E2 <- function(x, kap = a, rho.intvl = c(-0.8, 0.5)){
  re <- range.set(x1=x, rho.intvl = rho.intvl)
  cen <- (re[2]+re[1])/2
  amb <- re[2]-re[1]
  kap*cen+(1-kap)*amb
}
  #A1 <- A.f(sd.mat0, rho.vec, par.weight = a, i=1)
  cat("j = ",j, "\n")
for (i in seq_along(t.seq)){
  t <- t.seq[i]
  obj <- function(x) E2(x) - t*sum(mu.vec*x)
  opt.re <- optimize(obj, interval = c(-50,50))
  x1.opt <- opt.re$minimum
  x.opt <- c(x1.opt, 1-x1.opt)
  H.opt.mat[i, j] <- E1(x1.opt)
  mu.opt.mat[i, j] <- sum(mu.vec*x.opt)
}
}
## j =  1 
## j =  2 
## j =  3 
## j =  4 
## j =  5 
## j =  6 
## j =  7 
## j =  8 
## j =  9
#plot the matrix
#get the largest range
H.lim <- c(min(H.opt.mat), max(H.opt.mat))
mu.lim <- c(min(mu.opt.mat), max(mu.opt.mat))

for (j in seq_along(a.seq)){
  if (j==1){
    plot(H.opt.mat[,j],mu.opt.mat[,j], col=j,lty=j,
         type = "l", xlab = "H.opt", ylab = "mu.opt", 
         xlim = H.lim, ylim = mu.lim, 
         main = "Changing tuning par for Q1")
  } else {
    lines(H.opt.mat[,j],mu.opt.mat[,j], col=j,lty=j,
          type = "l")
  }
}

4.5 How does the optimal port depend on the paramter?

4.5.1 A1

4.5.1.1 sig1.intvl

#x1 vs sig1.intvl
x.seq <- sig1.low.seq <- seq(0.5, 4, .1)
y.seq <- sig1.amb.seq <- seq(0.1, 2, .1)

S <- function(a,b){
#both report the variance and current mean
  E1 <- function(x, lam = 0.5, rho.intvl = c(0.5, 0.5), 
                  sig1.intvl = lowamb2intvl(a,b), 
                  sig2.intvl = c(2, 2)){
  re <- range.set(x1=x, rho.intvl = rho.intvl, 
                  sig1.intvl = sig1.intvl, 
                  sig2.intvl = sig2.intvl)
  lam*re[2]+(1-lam)*re[1]
 }

 opt.re <- optimize(E1, interval = c(-50,50))
 x1.opt <- opt.re$minimum
 return(x1.opt)
}

z <- matrix(NA, nrow = length(x.seq), ncol = length(y.seq))

for (i in seq_along(x.seq)){
  for (j in seq_along(y.seq)){
    z[i,j] <- S(x.seq[i], y.seq[j])
  }
}

#z <- outer(x.seq, y.seq, S)

#fig1 <- plot3d(x.seq, y.seq, z.mat)
#rglwidget(elementId = "plot3drgl")

x <- x.seq
y <- y.seq
fig <- plot_ly(
  type = 'surface',
  # contours = list(
  #   z = list(show = TRUE, start = 0, end = max(z), size = 0.04)),
  x = ~x,
  y = ~y,
  z = ~z)
fig <- fig %>% layout(
    scene = list(
      xaxis = list(nticks = 10, title = "sig1.low"),
      yaxis = list(nticks = 10, title = "sig1.amb"),
      zaxis = list(nticks = 10, title = "x1.opt"),
      camera = list(eye = list(x = -1, y = -1, z = 0)),
      aspectratio = list(x = .9, y = .8, z = 0.6)))
fig

As \(\underline{\sigma}_1\) increases or the ambiguity of variance (of asset one) increases, the proportion on \(x1\) will decrease.

#x1 vs sig1.intvl
x.seq <- sig1.low.seq <- seq(0.5, 4, .1)
y.seq <- sig1.amb.seq <- seq(0.1, 2, .1)

S <- function(a,b){
#both report the variance and current mean
  E1 <- function(x, lam = 0.5, rho.intvl = c(-0.5, -0.5), 
                  sig1.intvl = lowamb2intvl(a,b), 
                  sig2.intvl = c(2, 2)){
  re <- range.set(x1=x, rho.intvl = rho.intvl, 
                  sig1.intvl = sig1.intvl, 
                  sig2.intvl = sig2.intvl)
  lam*re[2]+(1-lam)*re[1]
 }

 opt.re <- optimize(E1, interval = c(-50,50))
 x1.opt <- opt.re$minimum
 return(x1.opt)
}

z <- matrix(NA, nrow = length(x.seq), ncol = length(y.seq))

for (i in seq_along(x.seq)){
  for (j in seq_along(y.seq)){
    z[i,j] <- S(x.seq[i], y.seq[j])
  }
}

#z <- outer(x.seq, y.seq, S)

#fig1 <- plot3d(x.seq, y.seq, z.mat)
#rglwidget(elementId = "plot3drgl")

x <- x.seq
y <- y.seq
fig <- plot_ly(
  type = 'surface',
  # contours = list(
  #   z = list(show = TRUE, start = 0, end = max(z), size = 0.04)),
  x = ~x,
  y = ~y,
  z = ~z)
fig <- fig %>% layout(
    scene = list(
      xaxis = list(nticks = 10, title = "sig1.low"),
      yaxis = list(nticks = 10, title = "sig1.amb"),
      zaxis = list(nticks = 10, title = "x1.opt"),
      camera = list(eye = list(x = -1, y = 0, z = 0)),
      aspectratio = list(x = .9, y = .8, z = .6)))
fig
#x1 vs sig1.intvl
x.seq <- sig1.low.seq <- seq(0.5, 4, .1)
y.seq <- sig1.amb.seq <- seq(0.1, 2, .1)

S <- function(a,b){
#both report the variance and current mean
  E1 <- function(x, lam = 0.5, rho.intvl = c(-0.5, 0.5), 
                  sig1.intvl = lowamb2intvl(a,b), 
                  sig2.intvl = c(2, 2)){
  re <- range.set(x1=x, rho.intvl = rho.intvl, 
                  sig1.intvl = sig1.intvl, 
                  sig2.intvl = sig2.intvl)
  lam*re[2]+(1-lam)*re[1]
 }

 opt.re <- optimize(E1, interval = c(-50,50))
 x1.opt <- opt.re$minimum
 return(x1.opt)
}

z <- matrix(NA, nrow = length(x.seq), ncol = length(y.seq))

for (i in seq_along(x.seq)){
  for (j in seq_along(y.seq)){
    z[i,j] <- S(x.seq[i], y.seq[j])
  }
}

#z <- outer(x.seq, y.seq, S)

#fig1 <- plot3d(x.seq, y.seq, z.mat)
#rglwidget(elementId = "plot3drgl")

x <- x.seq
y <- y.seq
fig <- plot_ly(
  type = 'surface',
  # contours = list(
  #   z = list(show = TRUE, start = 0, end = max(z), size = 0.04)),
  x = ~x,
  y = ~y,
  z = ~z)
fig <- fig %>% layout(
    scene = list(
      xaxis = list(nticks = 10, title = "sig1.low"),
      yaxis = list(nticks = 10, title = "sig1.amb"),
      zaxis = list(nticks = 10, title = "x1.opt"),
      camera = list(eye = list(x = -1, y = 0, z = 0)),
      aspectratio = list(x = .9, y = .8, z = .6)))
fig

4.5.1.2 rho.intvl

#x1 vs sig1.intvl
x.seq <- rho.l.seq <- seq(-0.9, 0.9, .01)
y.seq <- rho.r.seq <- seq(-0.9, 0.9, .01)
load("z_rho.Rdata")
S <- function(a,b){
  if (a>b) {
    return(NA)
  } else {
    E1 <- function(x, lam = 0.5, rho.intvl = c(a,b), 
                  sig1.intvl = c(0.5,1), 
                  sig2.intvl = c(1,2)){
    re <- range.set(x1=x, rho.intvl = rho.intvl, 
                  sig1.intvl = sig1.intvl, 
                  sig2.intvl = sig2.intvl)
    lam*re[2]+(1-lam)*re[1]
 }
   opt.re <- optimize(E1, interval = c(-50,50))
   x1.opt <- opt.re$minimum
   return(x1.opt) 
    }
}

z <- matrix(NA, ncol = length(x.seq), nrow = length(y.seq))

for (j in seq_along(x.seq)){
  for (i in seq_along(y.seq)){
    z[i,j] <- S(x.seq[j], y.seq[i])
  }
}

#z <- outer(x.seq, y.seq, S)

#save(z, file = "z_rho.Rdata")
x <- x.seq
y <- y.seq
fig <- plot_ly(
  type = 'surface',
  # contours = list(
  #   z = list(show = TRUE, start = 0, end = max(z), size = 0.02)),
  x = ~x,
  y = ~y,
  z = ~z
)
fig <- fig %>% layout(
    scene = list(
      xaxis = list(nticks = 10, title = "rho.low"),
      yaxis = list(nticks = 10, title = "rho.up"),
      zaxis = list(nticks = 10, title = "x1.opt"),
      camera = list(eye = list(x = -1, y = 0, z = 0)),
      aspectratio = list(x = .9, y = .8, z = .6)))
fig

4.5.1.3 tuning parameter

(learn how to better interpret the results.)

5 Eigenvalue decomposition

(show the geometric result)

6 Directly visualize the function

6.1 Setup

#check sigma
d <- 1e-4
x.seq <- seq(d, 1-d,d)

# sd1 <- c(0.25, 1)
# sd2 <- c(0.75, 2)

sd1 <- c(0.5, 2)
sd2 <- c(2.5, 3)
sd1.c <- mean(sd1)
sd2.c <- mean(sd2)

s.set <- rep(sd1, each=2)/rep(sd2, 2)
rhol  <- -0.8

sd1[1]*sd2[2] < sd1[2]*sd2[1]
## [1] TRUE
rhol^2 > (sd1[1]*sd2[1])/(sd1[2]*sd2[2])
## [1] TRUE
rhol^2 < sd1[1]/sd1[2]
## [1] FALSE
# s.set.sort <- sort(c(s.set*(-rhol), 
#                      1/s.set*(-rhol)
#                      #s.set/(-rhol), 
#                      #1/(s.set*(-rhol))
#                      ))
# r.thres <- s.set.sort
# x.thres <- 1/(1+r.thres)

# r1 <- sd2[1]/(sd1[2]*(-rhol))
# r2 <- (sd1[1]*(-rhol))/sd2[2]
#important note
#the meaning of lower, upper in optim() is important!
re.mat <- matrix(NA, nrow = length(x.seq), ncol = 2)
val.seq <- numeric(length(x.seq))
for (i in seq_along(x.seq)){
 x1 <- x.seq[i]
 re <- optim(c(sd1.c,sd2.c), function(s) H(s[1],s[2], rho = rhol, x1 = x1), 
            lower = c(sd1[1], sd2[1]), 
            upper = c(sd1[2], sd2[2]), 
            method = "L-BFGS-B")
 #not lower=sd1, upper=sd2
 re.mat[i,]<- re$par 
 #get the min value
 val.seq[i] <- re$value
}
#any(is.nan(re.mat[,2]))
#re.mat
summary(re.mat[,1])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.5000  0.6668  2.0000  1.4400  2.0000  2.0000
summary(re.mat[,2])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   2.500   2.500   2.500   2.564   2.500   3.000
# optim(c(1,1,1), function(x) sum(x^2), lower = c(0,0,0), upper = c(2,2,2), method = "L-BFGS-B")
#r1 <- 1/abs(rhol)
#x.thres1 <- c(x.thres, 1/(1+r1))
matplot(x.seq, re.mat, type = "l", main = paste0("min cov ","\n", "with sd1 = ", "[", sd1[1], ",", sd1[2],"], ",
        "sd2 = ", "[", sd2[1], ",", sd2[2],"]"), lty = c(1,1))

abline(v = x.thres.f(s.set*(-rhol))[c(1,2)], lty = 2, col = "purple")

abline(v = x.thres.f(s.set/(-rhol))[c(1,3)], lty = 2, col = "blue")

# abline(v = x.thres.f((-rhol)/s.set), lty = 2, col = "brown")
# abline(v = x.thres.f(1/(s.set*(-rhol))), lty = 2, col = "pink")
#abline(v = x.thres.f(s.set/(-rhol)), lty = 1:4, col = 1:4)

#I have got it!
plot(x.seq, val.seq, type = "l")
abline(v = x.thres.f(s.set*(-rhol))[c(1,2)], lty = 2, col = "purple")
abline(v = x.thres.f(s.set/(-rhol))[c(1,3)], lty = 2, col = "blue")

t.break <- c(x.thres.f(s.set/(-rhol))[c(3,1)], 
             x.thres.f(s.set*(-rhol))[c(1,2)])

re.mat.c <- matrix(NA, ncol=2, nrow = length(x.seq))

sdl1 <- sd1[1]
sdr1 <- sd1[2]
sdl2 <- sd2[1]
sdr2 <- sd2[2]

q <- 1/(-rhol)
#check the function 
for (i in seq_along(x.seq)){
 x1 <- x.seq[i]
 x2 <- 1-x1
 r <- x2/x1 
 if (x1 <= t.break[1]){
   re.mat.c[i,] <- c(sdr1, sdl2)
 } else if (x1 <= t.break[2]) {
   re.mat.c[i,] <- c(r*sdl2/q, sdl2)
 } else if (x1 <= t.break[3]) {
   re.mat.c[i,] <- c(sdl1, sdl2)
 } else if (x1 <= t.break[4]) {
   re.mat.c[i,] <- c(sdl1, sdl1/(r*q))
 } else {
   re.mat.c[i,] <- c(sdl1, sdr2)
 }
 
#as x1 increases, r decreases 
 # if (r <= r.break[1]){
 #   re.mat.c[i,] <- c(sdl1, sdr2)
 # } else if (r <= r.break[2]) {
 #   re.mat.c[i,] <- c(sdl1, q*sdl1/r)
 # } else if (r <= r.break[3]) {
 #   re.mat.c[i,] <- c(sdl1, sdl2)
 # } else if (r <= r.break[4]) {
 #   re.mat.c[i,] <- c(r*sdl2/q, sdl2)
 # } else {
 #   re.mat.c[i,] <- c(sdr1, sdl2)
 # }
 
}


matplot(x.seq, re.mat.c, type = "l", main = paste0("min cov.check ","\n", "with sd1 = ", "[", sd1[1], ",", sd1[2],"], ",
        "sd2 = ", "[", sd2[1], ",", sd2[2],"]"), lty = c(1,1), ylim = c(0.5, 3))

abline(v = x.thres.f(s.set*(-rhol))[c(1,2)], lty = 2, col = "purple")

abline(v = x.thres.f(s.set/(-rhol))[c(1,3)], lty = 2, col = "blue")

r.seq <- (1-x.seq)/x.seq
matplot(log(r.seq), re.mat, type = "l", lty = c(1,1))

re.mat.d <- apply(re.mat, 2, diff)
matplot(x.seq[-1], re.mat.d, type = "l")

ind1 <- which(re.mat.d[,1]<0)
x.t1 <- c(min(x.seq[ind1]), max(x.seq[ind1]))
ind2 <- which(re.mat.d[,2]>0)
x.t2 <- c(min(x.seq[ind2]), max(x.seq[ind2]))
x.t <- c(x.t1, x.t2)
x.t
## [1] 0.5000 0.7999 0.8620 0.8823
#x.t*(-rhol)
r.t <- (1-x.t)/x.t
r.t
## [1] 1.0000000 0.2501563 0.1600928 0.1334013
r.t*(-rhol)
## [1] 0.8000000 0.2001250 0.1280742 0.1067211
r.t/(-rhol)
## [1] 1.2500000 0.3126953 0.2001160 0.1667517
# x.seq <- seq(.01, .99,.001)
re.mat <- matrix(NA, nrow = length(x.seq), ncol = 2)
for (i in seq_along(x.seq)){
 x1 <- x.seq[i]
 re <- optim(c(1,2), function(s) -H(s[1],s[2], rho = rhol, x1 = x1), 
            lower = c(sd1[1], sd2[1]), 
            upper = c(sd1[2], sd2[2]), 
            method = "L-BFGS-B")
 re.mat[i,]<- re$par 
}

#re.mat

matplot(x.seq, re.mat, type = "l", main = "max cov", lty = c(1,1))

#r.seq <- (1-x.seq)/x.seq
#matplot(r.seq, re.mat, type = "l")
#re.mat.d <- apply(re.mat, 2, diff)
#matplot(x.seq[-1], re.mat.d, type = "l")
tan.theta <- function(x1, rho = -0.5){
  x2 <- 1-x1
  a <- 2*x1*x2*rho
  D <- sqrt((x1^2-x2^2)^2 + a^2)
  b <- x1^2-x2^2 + D
  a/b
}
plot(tan.theta, from = 0, to = 1)

#tan.theta(0)
tan.theta(1)
## [1] 0
theta.f <- function(x1, rho = -0.5){
  if(x1==0) return(pi/2)
  if(x1==1) return(0)
  x2 <- 1-x1
  a <- 2*x1*x2*rho
  D <- sqrt((x1^2-x2^2)^2 + a^2)
  b <- x1^2-x2^2 + D
  atan(a/b)+pi
}
x <- seq(0,1,.01)
#sapply(x, theta.f)
plot(x, sapply(x, theta.f), type = "l")

#plot(theta.f, from = 0, to = 1)

theta.f(0.99999)
## [1] 3.141588
# px1 = 0.999
# x2 <- 1-x1
# a <- 2*x1*x2*rho
# D <- sqrt((x1^2-x2^2)^2 + a^2)
# b <- x1^2-x2^2 + D
# atan(a/b) 

6.2 Cen-amb graph

Center-ambiguity graph (to see the center and ambiguity of variance interval for a portfolio under different weights)

#check sigma
d <- 1e-4
x.seq <- seq(d, 1-d,d)

sd1 <- c(0.25, 1.5)
sd2 <- c(1, 2)

# sd1 <- c(0.5, 2)
# sd2 <- c(2.5, 3)

sd1.c <- mean(sd1)
sd2.c <- mean(sd2)

s.set <- rep(sd1, each=2)/rep(sd2, 2)
rhol  <- -0.8
rhor <- -0.2

sd1[1]*sd2[2] < sd1[2]*sd2[1]
## [1] TRUE
rhol^2 > (sd1[1]*sd2[1])/(sd1[2]*sd2[2])
## [1] TRUE
rhol^2 < sd1[1]/sd1[2]
## [1] FALSE
title.main <- paste0("sd1 = ", "[", sd1[1], ",", sd1[2],"], ", "sd2 = ", "[", sd2[1], ",", sd2[2],"]", ", ", 
                     "rho = ", "[", rhol, ",", rhor,"]")
print(title.main)
## [1] "sd1 = [0.25,1.5], sd2 = [1,2], rho = [-0.8,-0.2]"
#important note
#the meaning of lower, upper in optim() is important!
# re.mat <- matrix(NA, nrow = length(x.seq), ncol = 2)
val.min.seq <- val.max.seq <- numeric(length(x.seq))
for (i in seq_along(x.seq)){
 x1 <- x.seq[i]
 #min var
 re.min <- optim(c(sd1.c,sd2.c), function(s) H(s[1],s[2], rho = rhol, x1 = x1), 
            lower = c(sd1[1], sd2[1]), 
            upper = c(sd1[2], sd2[2]), 
            method = "L-BFGS-B")
 val.min.seq[i] <- re.min$value
 #max var
 re.max <- optim(c(sd1.c,sd2.c), function(s) -H(s[1],s[2], rho = rhor, x1 = x1), 
            lower = c(sd1[1], sd2[1]), 
            upper = c(sd1[2], sd2[2]), 
            method = "L-BFGS-B")
 val.max.seq[i] <- -re.max$value
}

matplot(x.seq, cbind(val.min.seq,val.max.seq), type = "l", lty = c(1,1), main = title.main)

#center
var.cen.seq <- (val.min.seq+val.max.seq)/2

plot(x.seq, var.cen.seq, type = "l",
        main = title.main)

#amb
var.amb.seq <- val.max.seq - val.min.seq

plot(x.seq, var.amb.seq, type = "l",
        main = title.main)

#plot it 
plot(var.cen.seq, var.amb.seq,
        main = title.main, type = "p")

#minimizer is very close to each other
matplot(x.seq, cbind(var.cen.seq,var.amb.seq), type = "l", lty = c(1,1), main = title.main)

abline(v = x.seq[which.min(var.cen.seq)], col = 1, lty=2)
abline(v = x.seq[which.min(var.amb.seq)], col = 2, lty=2)

library(ggplot2)
x1 <- x.seq
dat1 <- data.frame(x=var.cen.seq, y=var.amb.seq)
ggplot(dat1, aes(x=x, y=y)) + geom_point(aes(colour = x1)) + xlab("Var.center") + ylab("Var.amb") + ggtitle(title.main) +
scale_colour_gradient(low = "#00FFFF", high = "#000099")

#scale_colour_gradient(low = "lightblue", high = "darkblue")
#use different color to show the value of x1
#similar to efficient frontier 

It seems that there is a unique minimizer for smallest cen and ambiguity for a large class of objective function (linear combinations of center and ambiguity, use sum of square of two components). The minimizer does not depend much on the format of the function.

The mean part will also change monotonically w.r.t \(x1\), because \[ \mu = x_1\mu_1 + x_2\mu_2 \]

#write a function 
#plot.ambcen
plot.ambcen.var(sd1 = c(0.25, 1.5), sd2 = c(1,2), rho = c(-0.8, -0.2))

plot.ambcen.var(sd1 = c(0.25, 1.5), sd2 = c(1,2), rho = c(-0.8, 0.2))

plot.ambcen.var(sd1 = c(0.5, 2.26), sd2 = c(1,2.03), rho = c(-0.8, -0.4))

#we can notice some "discontinuity" 
#it should only come from some numerical issues, and the current x1 may go through some threshold. 
# plot.ambcen.var(sd1 = c(0.2, 1), sd2 = c(1.2,2.4), 
#                 rho = c(-0.8, -0.2))

(Do we also want to consider the mean part in this picture?)

#create a shiny app if needed 
#We can also change how the amb-cen curve changes depending on the parameter setup
#to let user change the parameter in a convenient way

#show the optimal x1 by drawing a straight line or using a circle (to better show the objective function)

library(shiny)

ui <- basicPage(
  plotOutput("plot1", click = "plot_click"),
  verbatimTextOutput("info")
)

server <- function(input, output) {
  output$plot1 <- renderPlot({
    plot(mtcars$wt, mtcars$mpg)
  })

  output$info <- renderText({
    paste0("x=", input$plot_click$x, "\ny=", input$plot_click$y)
  })
}

shinyApp(ui, server)

6.3 Min-max var plot

x1 <- x.seq
dat1 <- data.frame(x=val.min.seq, y=val.max.seq)
ggplot(dat1, aes(x=x, y=y)) + geom_point(aes(colour = x1)) + xlab("Var.min") + ylab("Var.max") + ggtitle(title.main) +
scale_colour_gradient(low = "#00FFFF", high = "#000099")

#it becomes harder to see
plot.ambcen.var(sd1 = c(0.25, 1.5), sd2 = c(1,2), rho = c(-0.8, -0.2),
                plot.minmax.var = TRUE)

plot.ambcen.var(sd1 = c(0.25, 1.5), sd2 = c(1,2), rho = c(-0.8, 0.2), 
                plot.minmax.var = TRUE)

re <- plot.ambcen.var(sd1 = c(0.2, 2), sd2 = c(1,1.8), rho = c(-0.5, 0.5),
                weight.obj = 0.5, plot.minmax.var = FALSE, 
                return.ind = TRUE, plot.ind = TRUE)

#check the break point
weight.seq <- seq(0.1, 0.9, .1)
x.opt.seq <- numeric(length(weight.seq))
for (i in seq_along(weight.seq)){
  w <- weight.seq[i]
  re <- plot.ambcen.var(sd1 = c(0.2, 2), sd2 = c(1,1.8), 
                        rho = c(-0.5, 0.5),
                weight.obj = w, plot.minmax.var = FALSE, 
                return.ind = TRUE, plot.ind = FALSE)
  x.opt.seq[i] <- re$x1.opt
}
plot(weight.seq, x.opt.seq, type = "l") #it will go through a dramatic change

re <- plot.ambcen.var(sd1 = c(0.2, 2), sd2 = c(1,1.8), rho = c(-0.5, -0.4),
                weight.obj = 0.5, plot.minmax.var = FALSE, 
                return.ind = TRUE, plot.ind = TRUE)

any(is.na(re$var.cen.seq))
## [1] FALSE
any(is.na(re$var.amb.seq))
## [1] FALSE
x.opt.seq2 <- numeric(length(weight.seq))
for (i in seq_along(weight.seq)){
  w <- weight.seq[i]
  re <- plot.ambcen.var(sd1 = c(0.5, 3), sd2 = c(1,3), rho = c(-0.5, -0.4),
                weight.obj = w, plot.minmax.var = FALSE, 
                return.ind = TRUE, plot.ind = FALSE)
  x.opt.seq2[i] <- re$x1.opt
}
plot(weight.seq, x.opt.seq2, type = "l") #it is not a dramatic change 

#the optimal portfolio is robust to the choice of the weight in the objective function (the sharper or narrower the nose is, the less sensitive it will become.)

6.4 Dynamic programming

#write a function by ourselves
h1 <- function(sig2, x1 = 0.5, rhol = -0.5, 
               sd1 = c(0.5, 2.0)){
  re <- optimize(function(sig1) H(sig1=sig1, sig2=sig2, 
                                  rho = rhol, x1 = x1), 
                 interval = sd1)
  re$objective
}

h1(2)
## [1] 0.75
plot(Vectorize(h1), from = 0.1, to = 5)

#create a shiny app if needed